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Abstract 



We present the results of first-principles molecular-dynamics simulations of 
molten silicates, based on the density functional formalism. In particular, the 
structural properties of a calcium aluminosilicate [ CaO-Al203-Si02 ] melt are 
compared to those of a silica melt. The local structures of the two melts are 
in good agreement with the experimental understanding of these systems. In 
the calcium aluminosilicate melt, the number of non-bridging oxygens found 
is in excess of the number obtained from a simple stoichiometric prediction. In 
addition, the aluminum avoidance principle, which states that links between 
AIO4 tetrahedra are absent or rare, is found to be violated. Defects such 
as 2-fold rings and 5-fold coordinated silicon atoms are found in comparable 
proportions in both liquids. However, in the calcium aluminosilicate melt, 
a larger proportion of oxygen atoms are 3-fold coordinated. In addition, 5- 
fold coordinated aluminum atoms are observed. Finally evidence of creation 
and anihilation of non-bridging oxygens is observed, with these oxygens being 
mostly connected to Si tetrahedra. 

PACS numbers: 61.20.Ja,71.15.Pd,91.60.-x,61.43.Fs 
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1. INTRODUCTION 



Silicate melts are the precursors to industrially and technologically important materials 
including ceramics and nuclear and industrial waste confinement glasses. They also occur 
naturally in the form of geologic magmas. Despite their technological relevance and geo- 
physical importance, however, their microscopic characteristics are not well understood [|TJ. 
The primary reason for this is that the structure of the melt is far more difficult to char- 
acterize than that of crystalline silicates, necessitating a combination of indirect methods. 
Most of the information about the liquid structures come from extrapolations of studies 
of glasses |2}^TB|, but large changes in properties, such as heat capacity, thermal expan- 
sion, compressibility, etc., as a result of heating suggest significant structural changes with 
increasing temperature [|T^-|TB[. 

In recent years, theoretical studies based on classical molecular-dynamics (MD) simu- 
lations have been able to treat binary or ternary silicate systems with reasonable success. 
However, accurate and reliable descriptions of systems containing more than three different 
atomic species by classical MD methods has proved far more difficult, although several stud- 
ies have been able to predict general trends that are consistent with experimental results 
T7H21~1. This problem is particularly acute when cations such as Na + or Ca 2+ are intro- 



duced in silicate glasses, as it is generally difficult to find interatomic potentials that can 
accurately treat both the covalent and the ionic nature of the interactions and are also ca- 
pable of describing the bond breaking and forming events that can occur in such chemically 
diverse environments. These difficulties can be circumvented by employing the ab initio 
molecular dynamics method, in which internuclear interactions are calculated "on the fly" 
from electronic structure calculations. Indeed, recent ab initio MD studies of silica glass and 
of its melt have demonstrated the ability of this approach to describe the local structure and 
dynamics of such systems with reasonable accuracy. These studies have also highlighted the 
advantages of employing a method which, additionally, allows direct access to the electronic 
properties of the system p2| , f23|l . 

It is well known that many important macroscopic properties of silicate melts, such as 
the viscosity, the glass transition temperature T g , or the resistance to chemical change, for 
example via corrosion, are dramatically altered by changes in the composition [p], p^ , p!8| . For 
instance, introduction of Na + ions into a SiC>2 melt causes the viscosity at a given T g /T ratio 
to decrease by several orders of magnitude from that of the pure SiC>2 [[|]. The presence of 
cations such as Na + , K + , Ca 2+ or Mg 2+ , known as network modifier cations, induces such 
changes by breaking some fraction of the Si-0 bonds thereby creating non-bridging oxygens 
and disrupting the tetrahedral silicate network. Non-bridging oxygens (NBO) are oxygen 
atoms which do not connect two tetrahedral cations or network-forming atoms, such as 
Si. Non-bridging oxygens provide relatively weak connections between the network forming 
atoms and the network modifier cations. However, when other network-forming atoms such 
as aluminum are introduced into the system, there is a gradual conversion of non-bridging 
oxygens into bridging oxygens. This arises from the fact that most of the Al atoms are 
tetrahedrally coordinated (AlOj), and the resulting negative charge compensates the pos- 
itive network modifier cation charge. In such cases, non-bridging oxygens can be created 
and the network broken only if there is an excess of network modifier cations, and it is for 
this reason that the viscosity of ternary liquids, such as CaSi0 3 or Na 2 Si0 3 , progressively 
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increases as the network modifier oxide (CaO or Na20) is replaced by AI2O3. The conven- 
tional explanation for the increase in viscosity is the transformation of non-bridging oxygens 
into bridging oxygens as the concentration of AI2O3 is increased. Generally, the number of 
NBO can be predicted based on a knowledge of the composition by simple stoichiometric 
arguments. However, it has recently been shown that such simple stoichiometric predictions 
are not exactly fulfilled in calcium aluminosilicate (CAS) glasses and that a small propor- 
tion of NBO can be present even if all the modifier cations should, in principle, exactly 



compensate the A10 4 tetrahedra fl4|T3 



In this paper, the results of an ab initio molecular-dynamics simulation of a calcium 
aluminosilicate (CAS) melt are presented and its microscopic characteristics are compared 
to those of a pure silica melt. To our knowledge, these are the first fully ab initio MD 
studies of the CAS melt. We have chosen a CaO-Al 2 03-Si02 system with a composition 
as close as possible to the basic composition of the confinement matrix for the nuclear 
wastes. This system also possesses a local structure close to those of some rapid cooling 
magmas. The chosen composition contains more Ca 2+ ions than are needed to compensate 
the AIO4 tetrahedra, thus leading to the formation of non-bridging oxygens. By carrying 
out a comparative study of the CAS and pure silica melts, it is possible to describe the 
detailed modifications in the network due to the presence of Al and Ca 2+ in the system. 

This paper is organized as follows: In Sec.|2|, the ab initio methodology is briefly described 
and the details of the particular simulations performed here are given. In Sec.^|, main results 
of the comparative study are presented, including structural properties of the CAS and silica 
systems. These results are discussed in Sec|| and conclusions are given in Sec.[|. 



2. SIMULATION DETAILS 



Equilibrated configurations of the two liquids were generated by classical molecular dy- 
namics runs, the details of which are described in Sees, p A| and [2 B| below, and were 
subsequently used to initialize the ab initio MD simulations. The two systems were then 
equilibrated within a Car-Parrinello (CP) ab initio MD run |24| performed with the ab ini- 
tio MD code, CPMD |[25|| . In the ab initio MD simulations, the electronic structure was 
treated via the Kohn-Sham (KS) formulation of density functional theory p6[ within the 
local density approximation for the pure silica system and within the generalized gradient 
approximation for the CAS system employing the B-LYP functional |27|,p8| . The KS orbitals 



were expanded in a plane-wave basis at the T-point of the supercell up to an energy cutoff 
of 70 Ry for both systems. Core electrons were not treated explicitly but were replaced 
by atomic pseudopotentials of the Bachelet-Hamann-Schluter type for silicon []29| and the 
Troullier-Martins type for oxygen pOf . A Goedecker-type pseudopotential |H] was employed 



for aluminum, and a Goedecker-type semi-core pseudopotential was employed for calcium. 
The choices of the pseudopotentials, exchange and correlation functionals and plane-wave 
cutoff are justified by previous studies carried out on amorphous Si0 2 |23[] as well as total 
energy calculations carried out on small molecules (see Table 0). 
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A. Molten Si0 2 



The molten silica system contains 26 SiC>2 units in a cubic box of edge length 10.558 A, 
which corresponds to a mass density of 2.2 g-cm -3 . The density of the glass was chosen so 
that the configurations could later be used in quenching runs to generate glass structures. 
Although the density is, therefore, a little too high compared to the real liquid, it is not 
expected that this will significantly affect our findings, which are based on the comparison 
of network and disrupted network systems. The SiC>2 initial configuration was obtained by 
melting a 216 SiC>2 units /3-cristobalite crystal at 7000 K with classical molecular dynamics 
using the van Beest, Kramer and van Santen (BKS) potential and then cooling it to 

4200 K using the same potential. At this temperature, a cubic box of edge length 10.558 A 
and containing 26 SiC>2 was extracted from the 216 SiC>2 system and equilibrated during ~ 
35 ps. 

The classically equilibrated SiC>2 liquid configuration was further equilibrated within a 
6-ps ab initio MD run at 4200 K using a time step of 0.096 fs, then quenched to 3500 K at 
a quench rate of 3 10 15 K-s -1 with the same time step, and finally equilibrated at 3500 K 
for 6 ps using a time step of 0.108 fs. In order to achieve rapid equilibration and efficient 
canonical sampling of the system, a separate Nose-Hoover chain thermostat |37| was placed 
on each ionic degree of freedom (known as "massive" thermostatting and an additional 



Nose-Hoover chain thermostat was placed on the electronic degrees of freedom |37];|39[]. In 
all fictitious electronic "mass" parameter, /i (having units of energy x time 2 ) of 600 

a.u. was employed. 



B. Molten calcium aluminosilicate 

The CAS system contains 22 SiC>2, 4 AI2O3 and 7 CaO, which gives approximately 67 %, 
12 % and 21 % molar percentages of these units, respectively, and a total of 100 atoms. For 
this particular composition, there are 8 Al atoms which give rise to 8 AIO4 tetrahedra under 
the assumption that all Al atoms form tetrahedra. Four of the Ca 2+ ions then compensate 
the negative charges of the AlOj, leaving three Ca 2+ that can break the network and create, 
ideally, 6 non-bridging oxygens. The system is confined in a cubic box with an edge length 
of 11.3616 A, which corresponds to a mass density of 2.4 g-cm -3 . This density has been 
chosen by extrapolating to 2500 K the data obtained by Courtial and Dingwell for a system 
of close composition |14}| . For this case, classical MD simulations were also carried out on 



systems containing 100 and 5184 atoms with the same composition described above. By 
comparing the structural properties obtained from the classical MD simulations at the two 
system sizes, it was possible to estimate the finite-size effects on the ab initio MD data and 



to validate the choice of the system size for the ab initio simulations [40|. 

The initial configuration of the CAS system for the ab initio MD simulation was generated 
using a Born-Mayer-Huggins potential [1|J in a classical MD run to obtain a melt at 2000 
K. The CAS liquid was then heated to 2300 K with CPMD and further equilibrated for 2 ps 
using a time step of 0.12 fs and an electron mass parameter of 800 a.u. It was then heated 
again to 3000 K and equilibrated for 6.8 ps, during which stronger diffusion effects occurred 
than at the lower temperature. Again, rapid equilibration and efficient canonical sampling 
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was achieved by coupling a Nose-Hoover chain thermostat p7| , |39f to each ionic degree of 
freedom. 



3. PRESENTATION OF RESULTS 

In this section, the structural properties of the CAS and Si02 melts are presented in terms 
of network pair correlation functions, angle distributions, examination of Al-O-Al linkages, 
proportion of NBO, and cation pair correlation functions. A full discussion of these results 
is presented in Sec. f|. 



A. Network pair correlation functions 

In this subsection, the pair correlation functions (PCF) corresponding to the network- 
forming atoms (Si, O and Al for CAS) are presented for the CAS melt and compared to 
those of the silica melt when appropriate. 

The Si-0 PCF and the corresponding integrated coordination number are almost identi- 
cal for the two systems, which confirms that the basic tetrahedral unit is conserved between 
these two liquids (Fig. |](b)). Moreover, the similarity in CAS between the Al-0 (Fig. |l|(d)) 
and Si-0 PCFs is clear evidence of the fact that the Al atoms can substitute for the Si 
atoms at the center of the tetrahedra. The ability of Al to replace Si in the network has 
also been observed in experimental studies of aluminosilicate glasses [[T|-|T]| • The slight shift 



of the Al-0 peak to a higher r value is consistent with previously observed and calculated 
Al-0 and Si-0 bond lengths in aluminosilicates as well as with the larger covalent radius 
of Al compared to Si. 

Comparison of the Si-Si and 0-0 PCFs between the CAS and Si02 liquids also shows 
a slight shift of the first peak toward higher r values in the CAS case, and the plateau in 
the running Si-Si coordination number is considerably lower in the CAS system. In order 
to explain the shift in the Si-Si peak, we first note that in silica, the first Si-Si neighbors 
correspond to neighboring tetrahedra (i.e. the first Si-Si distances correspond to the Si-O-Si 
linkages). In the CAS system, however, some of the oxygens, in particular, the NBO, are 
connected to only one network forming atom (Si or Al) and one network modifier (Ca), thus 
forming Si-O-Ca or Al-O-Ca linkages. Therefore, some of the first-neighbor Si-Si distances 
are due to these more complex linkages. The reduction in the Si-Si coordination is due 
primarily to the fact that some of the Si atoms are replaced by Al atoms. As a result of 
the Al substitution for Si and the presence of NBOs, there are fewer direct Si-O-Si linkages 
in CAS, leading to an average coordination of 2.4 compared to 4 in the pure silica system. 
The shift in the position of the first 0-0 peak in the CAS system is simply a reflection of 
the fact that the Al-0 bond length is larger than that of Si-0. Thus, if Al substitutes for Si, 
maintaining both the regular tetrahedral coordination and the angles between neighboring 
tetrahedra as in silica, then the observed shift in the 0-0 peak is expected. 
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B. Angles distributions 



Evaluation of the angles distributions and coordination numbers presented in this sub- 
section was based on distance cutoffs determined from the first minimum of the PCFs (2.38 
A for Si-0 and 2.56 A for Al-O). 

In order to investigate the effect of Al substitution for Si on the tetrahedral angles, the 
distributions of the O-Al-0 and the O-Si-0 angles for both systems were computed and 
are shown in Fig. |2|. The figure shows that the O-Al-0 distribution in CAS is somewhat 
broader than the O-Si-0 distribution in both CAS and in silica. The mean values of these 
distributions, however, are very similar (108.8 ° ± 15 for O-Si-0 in CAS, 108.3 ° ± 17 for 
O-Si-0 in Si02 and 107.6 ° ± 21 for O-Al-O), suggesting that the angles between tetrahedra 
remain approximately unchanged upon Al substitution. Therefore, the primary cause of the 
shift in the first peak of the 0-0 PCF is the increased Al-0 bond length. 

The relative broadness of the O-Al-0 angle distribution compared to that of O-Si-0 
likely signifies an Al-0 coordination that is different from 4. In order to check this assertion, 
the average oxygen coordination of the Al and Si atoms was computed. The histogram of 
the number of oxygen neighbors of Al and Si is shown in Fig. [3] (left panel). From the fig- 
ure, it can be seen that, although Al substitutes for Si in the network, the fraction of 3- fold 
and 5-fold coordinated Al exceeds that of Si. Recent high-temperature NMR measurements 
of aluminosilicate melts |16| , ^T| directly revealed the presence of 6-fold coordinated Al and 
strongly suggested the possibility of 5-fold coordinated Al, although the latter were not 
directly observed in these experiments. However, MAS NMR experiments [0,0 have pro- 
vided clear evidence of the existence of small amounts of both 5- and 6-fold Al coordination 
in binary Al 2 03-Si02 and in ternary CaO-Al 2 03-Si02 glasses [52,13]. Indeed, the average 
proportion of 5-fold coordinated Al in the present simulations is also small compared to 
4-fold coordinated Al. However, we observe less than 1 % of 6-fold coordinated Al for this 
particular composition. 

We also investigated the effect of the Al substitution for Si on the angles between neigh- 
boring tetrahedra. The Si-O-Al and Si-O-Si angles distributions were computed for the 
CAS melt and compared to the Si-O-Si angle distribution for the Si0 2 melt (Fig. [|). The 
Si-O-Si distributions are very similar in the two systems, only showing a slight difference 
in the intensity of the shoulder around 90 °- 100 °. It can be shown that values of the 
Si-O-Si angles around 90 °- 100 ° can be attributed to 2-membered rings and/or to oxygen 
tri-clusters, oxygens bound to three network- forming atoms (see Sec. |j). A comparison of 
the proportions of these units in the two different systems is consistent with the differences 
observed in the angle distributions and will be discussed in more detail in section f|. In CAS, 
the Si-O-Al angle distribution shows a higher shoulder around 90 °- 100 ° than the Si-O-Si 
one. Again the Si-Al 2-membered rings and the oxygen tri-clusters are at the origin of this 
shoulder. Overall, the introduction of Al and Ca does not significantly affect the Si-O-Si 
angle distribution in the molten state. 



C. Al-O-Al linkages 

Although the Al-Al PCF (Fig. |l](e)) is somewhat noisy due to the small number of 
Al atoms in the system, this correlation function exhibits a first peak at ~ 3.2 A, which 
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indicates that some Al-O-Al linkages are present in the system. This would appear to 
be in direct violation of the so called Al avoidance principle or Lowcnstcin's rule [44|. an 



empirical rule which states that two Al tetrahdera are never found linked by an oxygen atom 
in aluminosilicate crystals and glasses, at low Al content. However this rule has been found 
experimentally not to be exactly fulfilled in some glasses and melts, in particular when Ca 
atoms are present, leading to Al/Si alternance disorder |45|-|47[|. 

The degree of Al avoidance violation can be quantified by examining the average number 
of Al-O-Al linkages formed over the course of the simulation. By computing the number of 
Al atoms around each oxygen atom (using the first minimum of the g(r), 2.56 A, as the Al-0 
cutoff distance) and counting one Al-O-Al linkage when two Al atoms are found, an average 
value of 2.26 Al-O-Al linkages is obtained, i.e., ~ 57 % of Al atoms form Al-O-Al linkages, 
if they do not form chains. Moreover, half of the oxygen atoms involved in the Al-O-Al 
linkages are found to be 3-fold coordinated, on average. The proportion of Al-O-Al linkages 
is greater than that which would be obtained from a purely random model P5| . Thus, the 
Al atoms appear to favor Al-O-Al linkages in the CAS liquid and adhere only minimally to 
the Al avoidance principle. That the presence of Ca 2+ cations preferentially favor Al-O-Al 
linkages and, thus, strong violation of the Al avoidance principle, has also been suggested 
by static ab initio calculations of clusters [ESI and is likely due to the greater aggregation 



of negative charge around the Al-O-Al linkage. In a disordered (liquid or glass) state, the 
relatively large Ca 2+ charge is, therefore, able to induce formation of such linkages. 



D. Non-bridging oxygens 

In the present CAS system, stoichiometry dictates that four of the Ca 2+ ions should 
compensate the eight aluminum tetrahedra (assuming all the Al are four- fold coordinated), 
leaving three Ca 2+ cations to create six non-bridging oxygens. The number of NBOs in the 
system can be computed by counting the number of oxygen atoms which have only one Si 
or Al neighbor and one Ca neighbor. The neighbors are defined to be any atoms within 
a sphere of radius determined by the first minimum of the corresponding PCF (2.38 A for 
O-Si, 2.56 A for O-Al, and 3.40 A for O-Ca - see Sec. and Sec. centered on 

each oxygen atom. The distribution of the number of non-bridging oxygens found in the 
CAS system during the simulation is depicted in Figure ||. The distribution is peaked at 
9, while the probability of the system's possessing six NBOs is relatively small. Thus, the 
average number of NBOs in the system is larger than would be predicted from a simple 
stoichiometric argument. In a recent experiment of Stebbins and Xu |12j| , it was found that 



in a particular CAS glass system in which the Ca 2+ ions perfectly compensate all of the 
negative Al tetrahdera, some NBOs are present, although stoichiometric arguments would 
predict that there should be none. Since the present study is concerned with the liquid state, 
where bond-length fluctuations and other dynamical effects are more significant than in the 
glass, we cannot directly compare our result with that of Stebbins and Xu. However, given 
that the trend is toward a larger number of NBOs in the glass than the stoichiometry would 
predict, the fact that our calculations predict a similar trend in the liquid accords well with 
the experimental result. 

In Figure H we present the first peaks of the PCFs between the oxygen atoms and the 
network- forming atoms (Si and Al), evaluated separately for bridging (BO) and non-bridging 
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(NBO) oxygen atoms. The maximum intensity of the X-BO and X-NBO (X=Si,Al) peaks 
are located at different r values, the X-NBO distances being shorter than the corresponding 
to those for X-BO. This result is in agreement with experimental results concerning sodium 
silicates and aluminosilicates Hf49| and with ah initio calculations on clusters [^8|,^0| . In the 
simulated CAS melt, the most probable distances are r(Si — BO) ~ 1.64 A and r(Si — NBO) 
« 1.58 A, r(Al - BO) « 1.75 A, and r(Al - NBO) « 1.70 A. It is also interesting to note 
that the Si-0 and Al-0 distances in the molten state are very close to their respective values 
in the glass. 



E. Calcium pair correlation functions 

In this subsection, structural properties involving the calcium atoms in the CAS melt, 
such as PCFs and coordination numbers, are presented. 

The calcium PCFs are depicted in Fig. [7]. The Ca-0 radial distribution function 
(Fig. |7|(b)) exhibits a first peak at approximately 2.33 A, in agreement with experimen- 
tal values obtained for glasses and minerals of similar composition The coordination 
of calcium atoms is found to be equal to 6.2 ± 1.3 on average (see Fig. §), while the ex- 
perimental value in calcium aluminosilicate glasses is estimated to be between 5 and 6 . 
This comparison is only meant as a qualitative one, since direct comparison between the 
liquid and glass states is not possible. Indeed, the high temperature of the molten state 
gives rise to large fluctuations in the Ca-0 coordination number. 

The Ca-Si and Ca-Al PCFs ((Fig.^(a) and (c)) also show well defined first peaks at ~ 
3.45 A, which are due either to direct bonds with non-bridging oxygens, i.e., Si-O-Ca and Al- 
O-Ca linkages, or to the proximity of Ca 2+ ions to compensate the negative AlOj groups. 
It can be shown that the Ca-Si peak is due to the former and the Ca-Al to the latter. 
Indeed, we observe mostly Si-NBO bonds and very few Al-NBO bonds (~ 91 % of the 
NBOs are connected to Si atoms), a result that is consistent with recent X-ray experiments 
on glasses of similar composition [[|. In Fig. [| in which PCFs are evaluated separately 
for bridging (BO) and non-bridging (NBO) oxygen atoms, we observe that Si PCFs possess 
similar characteristics for both BO and NBO. In contrast, the Al-NBO PCF is generally 
small compared to the Al-BO PCF, which indicates that almost all the NBO atoms are 
connected to Si tetrahedra. Thus, the first peak in the Ca-Al PCF is due to nearby Ca 2+ 
ions which, in this close proximity, are able to compensate the negative AlOj groups. Note 
that the Ca-Ca PCF, shown in Fig. |7|(d) does not exhibit a well defined first peak, which 
suggests that at liquid conditions, there is no ordering of the network modifier cations. 



4. DISCUSSION 



In order to explain the presence of excess non-bridging oxygens, Stebbins and Xu fl2 
proposed two possible structural units: AIO5 groups and tri-clusters. Tri-clusters are oxygen 
atoms bonded to three network-forming atoms, either Al or Si, and they can be of four types 
: oxygens bonded to three silicon atoms (3 Si), to two silicon and one aluminum atoms (2 Si 
- 1 Al), to one silicon and two aluminum atoms (1 Si - 2 Al), and to three aluminum atoms 
(3 Al). In the molten state we observed that the fraction of 3-fold coordinated oxygen atoms 



S 



is not negligible. A more detailed analysis of these tri-clusters in CAS showed that, due to 
the small number of aluminum atoms in the system, 3-A1 tri-clusters are absent and that 
the most numerous tri-clusters are the 2 Si - 1 Al and 2 Al - 1 Si. When added together and 
averaged over the trajectory, the percentage of oxygen atoms forming tri-clusters is around 
6.9 % in the CAS system, with large fluctuations of about ± 2.8 %. This result accords well 
with the idea of Stebbins and Xu |12| that a given number of oxygen tri-clusters should be 
present in the glass in order to compensate the formation of the non-bridging oxygens. It 
could be argued that the relatively small number of tri-clusters may be a consequence of 
thermal fluctuations at liquid conditions. Although thermal induced tri-cluster formation 
cannot be ruled out, it is worth noting that only 4.7 % of tri-clusters could be identified in 
the silica melt at higher temperature. 

On the other hand, in CAS, we also find a relatively high number of AIO5 units: on 
average, approximately ~ 1.9 or 23% of Al are 5-fold coordinated, however the fluctuations 
are such that in any configuration, anywhere between and 3 AIO5 units may exist. This 
result suggests that the presence of highly coordinated aluminum atoms could favor the 
creation of excess NBO atoms. Most experimental studies PHTTIj only show evidence of 4-fold 
coordinated aluminium atoms in aluminosilicate glasses. However, as discussed in Sec. p B| , 
experimental evidence of higher Al coordination numbers (mainly 5 and 6-fold) in the glass 
and molten state of calcium aluminosilicate systems exists, and, recently, evidence of high 



Al coordination in magnesium aluminosilicate glasses has been reported [04, 16 , 41-43 



In Fig. |5|, we present the distribution of the number of NBOs found during the simulation 
of the CAS melt. The relatively large width of the distribution indicates that the identity 
of the oxygen atoms (BO or NBO) does not remain constant during the simulation. Indeed, 
bond-breaking events that lead to the creation and the anihilation of non-bridging oxygens, 
i.e. Q4 and Q 3 exchanges, are observed. This mechanism has been suggested as underlying 
the induction of shear flow and, hence, a decrease in the viscosity in these systems |14]-|16|. 
Since almost no NBOs are found to be connected to Al atoms on average during the sim- 
ulation (see Sec. |3 FJ ), it is highly probable that more Al-0 bonds are broken than Si-0 
bonds. In silica, however, since there are no NBOs, a different mechanism must underly 
the flow process. Unfortunately, we can not extract direct dynamical quantities from the 
present simulations. This is mainly due to the fact that at high temperature in the molten 
state, the electronic gap is too small compared to k^T to ensure the decoupling of the ionic 
and the electronic degrees of freedom, which is needed for ab initio molecular dynamics in a 
constant energy ensemble. The use of thermostats becomes compulsory and the direct access 
to dynamical properties is no longer available. We are currently developing new techniques 
to treat this problem pHJ. 

Nevertheless, the structural characteristics of the melt can give some insight into its 
dynamical properties. For instance, it has been suggested that five-coordinated network 
cations could act as transition states for the flow process []T3|-[n|[43|] . This atomic-scale flow 



step was described in Ref. [ 13 1 as an oxygen atom changing from bridging to non-bridging, 



or vice versa: 'A non-bridging oxygen bonded to a modifier cation can approach a silicon 
atom and make it over-coordinated. Dissociation of the over-coordinated results in an 
exchange of the roles of oxygen from bridging to non-bridging In the present simulations, 
we find a significant fraction of 5-fold coordinated silicon atoms (~ 9.9 % in Si02 at 3500 K 
and ~ 3.8 % in CAS at 3000 K) and a large fraction of 5-fold coordinated aluminum atoms 
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in CAS (see Fig.|3|). These numbers show very large fluctuations around their average values 
which supports the idea that these units participate in the transition state of the flow steps. 
The relatively large number of 5-fold coordinated network cations in CAS compared to Si02, 
even at a lower temperature, could be a signature of a lower viscosity in the aluminosilicate 
melt. It is also interesting to note that the number of 5-fold coordinated Si atoms is larger 
in Si02 than in CAS. In the latter system, the 5- fold coordinated Al atoms, which are 
energetically more favorable, replace the 5-fold Si atoms as the transition state of the flow 
step. 

In Sec. |3], it was seen that the local order in the two silicate melts is only slightly affected 
by temperature and by the introduction of a small number of modifier cations. The Si 
tetrahedral units remain the most probable basic units of the system, with probabilities of 
85 % in Si02 and 89% CAS at 3000 K. From structural characteristics such as pair correlation 
functions or angle distributions, therefore, it is almost impossible to discern the disruption 
of the silicate network by the modifier cations, at high temperature. The effect can be seen, 
however, by looking at characteristic patterns in the network. An example is the oxygen- 
oxygen coordination number, which is a more sensitive probe of the network disruption. We 
have evaluated these coordination numbers by counting the number of oxygen neighbors 
around each oxygen atom inside a sphere of radius determined by the first minimum of 
the 0-0 radial distribution function. Histograms of the oxygen-oxygen coordinations are 
depicted in Fig. ^| for the two different systems. For Si02 at 3500 K, the distribution of 
the O coordination peaks around 7 or 8 atoms, which indicates that the oxygen atoms have 
more oxygen neighbors than would be expected from two connected tetrahedra. In the Si02 
glass, in which the network is not broken, the maximum of the distribution occurs at 6 which 
corresponds to the number of oxygen neighbors belonging to two connected tetrahedra. The 
fact that, at 3500 K, the maximum is located between 7 and 8 in silica shows that part of 
the network has been broken. One of our hypotheses is that 2-membered rings and highly 
coordinated network formers are the cause of the increased number of oxygen neighbors. 
Indeed, the percentage of 2-membered rings in silica is not negligible at 3500 K (~ 12.4 %). 
It could be argued that our liquid silica system is too compressed (its density is equal to 
the experimental density of the silica glass at room temperature) and that the high pressure 
induces the shift in the 0-0 correlation and the large number of edge-sharing tetrahedra. 
However, results of classical-molecular dynamics on liquid silica at several densities and 
temperatures show that this is not the case |]52| , |53l . At ~ 3200 K and zero pressure, the 0-0 
distribution ressembles the distribution in Fig. ^| for silica [52"] whereas at higher pressure 
(for a density of 2.35 g.cm~ 3 ) and temperatures up to 4200 K, the 0-0 distribution is 
peaked between 6 and 7, which is close to that of the glass structure [53fl . The shift of 
the 0-0 distribution with pressure is accompanied by a decreasing number of edge-sharing 
tetrahedra. This "non intuitive" behavior has already been seen in previous Monte Carlo 
simulations of liquid silica under pressure, in which it has been shown that the proportion of 



small rings decreases with increasing pressure |5J] . Given the similarities between our result 
and that of Ref. ||52|| , and given the small thermal expansion of silica (0.54 10 -6 K" 1 at 1000 



C |j55|)) it is clear that the high 0-0 correlation is a consequence of the high temperature 



of the system. 

The distribution of 0-0 coordination in CAS at 3000 K is shifted towards higher numbers 
compared to that of Si0 2 (Fig. 0). In the former, the disruption of the network is not only 
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due to the NBOs but also to the presence of highly coordinated Al atoms which weaken 



the structure | lq| . The CAS melt density is only slightly larger than that given by the 



extrapolation to 3000 K of experimental data |IJJ, therefore the high correlation of the 



oxygen atoms is clearly due only to the high temperature. The relative shifts of the 0-0 
distributions between the two different systems can be related to the disruption of the three 
dimensional network. In the two melts, thermal effects create defects such as 2-membered 
rings and SiOs units in comparable proportions. In CAS, however, the disruption of the 
network is also driven by the presence of non-bridging oxygens, which break some of the Si- 
O bonds, as well as by the presence of the aluminum atoms which possess high coordination 
numbers. Indeed, the presence of aluminum atoms has already been proven to be responsible 
for the increased fragility of the silicate network which can be traced to the fact that 



the Al have broader coordination distributions, leading, therefore, to an increased flexibility 
of the network structure than in Si02- 

The use of the Nose-Hoover chain thermostatting method |J7| in the present simulations 
allows the heat capacity at constant volume, C v to be computed efficiently using the relation: 



AE = Jk B T 2 C v (1) 



where AE = y (E 2 ) — (E) 2 denotes the fluctuations of the potential energy plus the kinetic 
energy of the ions over the course of the simulation. In silica at 3500 K, we find that the 
heat capacity at constant volume is equal to 91.5 J mol -1 K _1 . As a first approximation, 
we can roughly estimate the difference between C v and C p , the heat capacity at constant 
pressure, which is the experimentally measured quantity: 

C p -C v = TV j (2) 

where a is the thermal expansion coefficient and (5 the isothermal compressibility. Using the 
experimental values a=0.54 10" 6 K -1 , valid at 1000 ° C, and (3 = 2.74 10" 2 GPa" 1 at 300 
K we find a difference of C p - C v = 1.0 10" 3 J mol" 1 K _1 , so C p » 91.5 J mol" 1 K _1 . 
This value does not disagree with an extrapolation of the experimental values of C p from 
Ref. |56|] and is close to the values obtained by Scheidler et al. |57| from molecular dynamics 



simulation of the Si02 melt in the same temperature range, using the BKS potential [|35| 
Scheidler et al. computed the frequency dependent specific heat from T=2750 K up to 
T=6100 K and found a reasonable extrapolation of the experimental C P (T) above T g . 

For the CAS melt at 3000 K, we find a value of C v = 99.28 J mol -1 K _1 and an estimate 
of the difference, C p - C v = 0.189 J mol" 1 K~\ which gives C p w 99.47 J mol" 1 K~\ The 
C p - C v difference has been estimated using experimental values for calcium aluminosilicate 
systems of close compositions: a = 54.5 10~ 7 K _1 at 815 ° C (25 mol% CaO, 15 mol% 
A1 2 3 , 60 mol% Si0 2 ) H and (3 = 1.28 10~ 2 GPa" 1 at 300 K (26.7 mol% CaO, 13.3 mol% 



A1 2 3 , 60 mol% Si0 2 ) [0. Although this estimate of C p is relatively crude, it can be used to 
give an order of magnitude for C p at high temperature. To our knowledge, no experimental 
values of C p have been reported for a CAS system of close composition, at high temperature. 
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5. CONCLUSION 



The structural properties of two different molten silicates have been analyzed using first- 
principles molecular-dynamics. It has been observed that, even in the calcium aluminosili- 
cate melt, the basic structural unit of the system, i.e. the SiC>4 tetrahedron is not destroyed 
at high temperature and that the aluminum atoms can substitute for the silicon atoms in 
the center of the tetrahedra. Analysis of the structure of two melts show that the tempera- 
ture effects induce the creation of 2-membered rings, 3-fold coordinated oxygen atoms and 
5-fold cooridnated silicon and aluminum atoms. In CAS, in particular, a relatively large 
proportion of Al atoms are found to be 5-fold coordinated. 

In the CAS system, a larger number of non-bridging oxygens than would be predicted 
from stoichiometry was found. This excess of NBO is in agreement with recent experimental 
results on glasses of composition similar to the present composition [12] and is likely due to 
the presence of AIO5 units as well as oxygen tri-clusters (mainly 2 Si - 1 Al). Also in the 
CAS system, it is found that the aluminum avoidance principle is violated and that most of 
the non-bridging oxygens are located on silicon tetrahedra, both results being in agreement 
with recent experimental data P P5|J^6[] . Finally, evidence of the creations and anihilations 
of non-bridging oxygens has been obtained. These processes likely play a key role in the 
flow mechanism in this system. A detailed understanding of the diffusion process, although 
extremely interesting, would require considerable computer time and would represent a 
methodological challenge for the ab initio approach. While we are currently developing new 
techniques to treat this problem |5l] , the study of the dynamical properties in these systems 
remains, for the present, beyond the scope of this paper. In addition, since it is not possible 
to determine which of the dominant structural motifs in CAS are a direct consequence of 
thermal fluctuations, we intend, as future work, to perform a quench of the CAS system to 
room temperature in order to study its structural properties in the glassy state. 
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FIGURES 
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FIG. 1. Pair correlation functions of the network-forming atoms in the CAS system (bold 
lines) and in the silica system (thin lines). The integrated coordination numbers are also shown in 
dashed lines (bold and thin for CAS and Si02, respectively). 
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FIG. 2. Distribution of the O-Al-0 angles in CAS (dashed line) and of the O-Si-0 angles in 
CAS (bold line) and in Si0 2 (thin line). 
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FIG. 3. Left graph: Histograms of Pxo(N) (X = Si,Al), the probability of finding N oxygen 
atoms around the silicium (black) and aluminum (white) atoms in CAS. Right graph: Histogram 
of Psio(N), the probability of finding TV oxygen atoms around the silicium atoms in SiO"2. 
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FIG. 4. Distribution of the Si-O-Si angles in CAS (bold line) and in Si02 (thin line) and 
distribution of the Si-O-Al angles in CAS (dashed line). 
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FIG. 5. Average distribution of the number of NBO atoms in the CAS system. The size of the 
symbols indicate the maximum error bar associated with the choice of distance cutoff for Al-0 and 
Ca-O. 
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FIG. 6. Left graph: First peaks of the Si-BO (dashed line) and Si-NBO (bold line) pair corre- 
lation functions in CAS. Right graph: First peaks of the Al-BO (dashed line) and Al-NBO (bold 
line) pair correlation functions in CAS. 



17 



3.0 
2.0 
1.0 
0.0 
2.0 
1.0 
0.0 



(a) Si-Ca 




(c) Al-Ca 



(b) Ca-0 

h 


I 

/ 

1 

I I I 1 I I I 1 I I I 1 I I I 


(d) Ca-Ca 

... I ... I . . J L 


/ 

/ 

1 / 
Jm / 

V / 'V 

/ 

/ 



0.0 1.0 2.0 3.0 4.0 5.0 6.0 0.0 1.0 2.0 3.0 4.0 5.0 6.0 

r (A) r (A) 

FIG. 7. Pair correlation functions of the calcium atoms in the CAS system and the correspond- 
ing integrated coordination numbers. 
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FIG. 8. Probability distribution, Pc a o(N), for finding iV oxygen atoms around a calcium atom 
in the CAS system. The error bars were estimated by bracketing the first minimum of gca-o{ r ) 
between statistical lower and upper bounds and varying the cutoff distance between these bounds. 
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number of oxygen neighbors 

FIG. 9. Probability distribution, Poo(N), for finding N oxygens around a given oxygen atom 
in CAS and SiO"2. The error bars were estimated by bracketing the first minimum of go-o( r ) 
between statistical lower and upper bounds and varying the cutoff distance between these bounds. 
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TABLES 



TABLE I. Results of distances and angles obtained from total energy calculations carried out 
on AIH2OH, CaH 2 and CaH using the pseudo-potentials and the energy cutoffs given in the text. 
A Trouiller-Martins pseudo-potential was used for hydrogen. For these calculations, the recipro- 
cal-space cluster boundary conditions method of Martyna and Tuckerman was employed |32|] 



AlH 2 OH 


This work 


Ref. |33] 


Al-0 (A) 
Al-Hx (A) 
Al-H 2 (A) 

0-H3 (A) 

HiATh 2 

A10H3 


1.73(5) 
1.56(4) 
1.57(1) 
0.97(6) 
125.5 
122.3 


1.711 
1.565 
1.575 
0.956 
123.9 
122.3 


CaH 2 


This work 


Ref. 1 34] 


Ca-Hi (A) 
Ca-H 2 (A) 


2.03(0) 
2.03(1) 


2.00 
2.00 


CaH 


This work 


Ref. [34] 


Ca-H (A) 


1.97(5) 


2.00 
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